Taller INLA covid
1 Analisis Predictivo
El análisis predictivo emplea datos históricos para predecir eventos futuros. Normalmente, los datos históricos se utilizan para crear un modelo que capture diversos patrones temporales. Este modelo predictivo se usa entonces con los datos actuales para predecir lo que acontecera. Existen diversas metodologias para llevar acabo predicciones, en este taller optamos por el enfoque bayesiano para llevar acabo la prediccion.
1.1 Enfoque bayesiano : Distribucion predictiva aposteriori
El enfoque bayesiano, en el contexto del analisis predictivo, enfoca su analisis en el calculo de la distribucion predictiva aposteriori. Esta distribucion emplea la informacion de los datos asi como la distribucion aposteriori de los parametros de interes . Empleando estos insumos la distribucion predictiva podemos generar datos futuros; es decir, datos faltantes en un determinando horizonte.
\[\begin{align*} f(y_{pred}|y_{obs}) = \int_{\theta} f(y_{pred},\theta|y_{obs})d\theta =\int_{\theta} f(y_{pred}| \theta y_{obs})f(\theta|y_{obs})d\theta \end{align*}\]
1.2 INLA
Actualmente existen diversas metodologias para llevar acabo el calculo numerico de esta distribucion las cuales en R corresponden al empleo de un determinado paquete en. Algunos paquetes ampliamente conocidos son los siguientes :
- STAN
- Winbugs
- Jags
Relativamente hace no mucho tiempo la metolodogia de la Aproximación anidada integrada de Laplace o INLA(por sus siglas en ingles) ha cobrado importancia por su relativa eficiencia para estimar modelos complejos.
1.3 Prediccion como datos faltantes
El analisis de prediccion como en la practica puede ser entendida como predecir datos desconocidos en un determinado horizonte. En ese sentido estos datos a futuro pueden ser interpretados como datos faltantes o missings values, los cuales son el objeto de la predicción.
| Año | Semana | Numero de muertes |
|---|---|---|
| 2019 | 1 | 0 |
| 2019 | 2 | 1 |
| 2019 | 3 | 4 |
| 2020 | 1 | 0 |
| 2020 | 2 | 1 |
| 2020 | 3 | 4 |
| 2021 | 1 | NA |
| 2021 | 2 | NA |
| 2021 | 3 | NA |
Para llevar acabo predicciones empleando INLA, una aproximacion es justamente la anterior. Una vez dispuestos los datos como en la tabla anterior, basta con ajustar un modelo tipico en INLA, ya que INLA internamente calculará la distribucion predictiva aposteriori .
2 Analisis Descriptivo
library(tidyverse)
library(INLA)
library(yardstick)
library(gt)
library(innovar)
library(spdep)
library(reshape2)
db <- readRDS("db_excess_proc.rds")
db.lima <- db %>%
filter(prov=="LIMA")
db.lima %>%
ggplot()+
geom_line(aes(x=date,y=n),color="#011f4b",lwd=1.5)+
geom_point(aes(x=date,y=n),color="black",lwd=2.5)+
ylab("Numero de muertes")+
xlab("")+
ggtitle("Evolucion del Numero de muertes ")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5,size = "22"),
strip.background =element_rect(fill="#011f4b"),
strip.text = element_text(color = "white",face = "bold",size = "5pts")) 
2.1 Modelamiento con INLA
Para llevar acabo el ajuste de modelos empleamos la funcion inla(), la cual tiene diversos parametros, algunos de los mas importantes son :
data: Un objeto tipicamente de la clasedata.frameel cual representa los datos para ajustar el modeloformula: Un objeto de la clase formula que especifica la ecuacion que pretendemos ajustar como por ejemploinla(target ~ 1 + x1). En la formula podemos especificar efectos lineales (introduciendo el nombre de lavariable) o no lineales (empleandof()).verbose: Una variable del tipoboolean, la cual señala si se desea mostrar en consola el proceso de convergencia
lima.m.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.lima
)Los parámetros antes detallados son los esenciales para ejecutar el ajuste de un modelo empleando INLA.Sin embargo, algunos parametros extra a tener en consideracion son los siguientes :
family:Objeto de clasecharacter.Este parametro es crucial, pues determina la distribucion de la variable objetivo, por defecto se encuentra enfamily=Gaussian.control.compute:Objeto de la claselistpermite especificar el calculo de criterios de informacion tales comoaic,dic,waic.
lima.m.m1 <- inla(n ~ 1 + annus,
verbose = F,
data = db.lima,
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
)2.2 Prediccion: Enfoque Temporal
En el analisis de prediccion (forecast) tipicamente se emplea datos en el formato de series de tiempo para llevarla acabo. Para ello comunmente se asume que las series de tiempo pueden ser descompuestas del siguiente modo :
\[\begin{align*} y_{t} = S_{t}+T_{t}+e_{t} \end{align*}\]
donde :
\(S_{t}\) : Estacionalidad
\(T_{t}\) : Tendencia
\(e_{t}\) : Ruido blanco
2.2.1 Ajuste y prediccion
Los componentes antes detallados pueden ser modelados en el contexto de INLA empleando distribuciones apriori para cada componente tales como los siguientes:
- linear :
variable
- AR1 (Proceso autoregresivo de 1er orden) :
f(variable, model = "ar1")
- RW1 (Paseo aleatorio de 1er orden) :
f(variable, model = "rw1")
- RW2 (Paseo aleatorio de 2do orden) :
f(variable, model = "rw2")
Dado el comportamiento que exibe la series que prentedemos modelar, para los componentes :
- Tendencia : Emplearemos una distribucion apriori del tipo autoregresivo(AR1) o linear
- Estacionalidad : Emplearemos una distribucion apriori de tipo de paseo aleatorio de 1er o 2do orden
Tomando en cuenta estas breves sugerencias procedemos a ajustar modelos para ello en primer lugar ,tomando en cuenta lo expuesto anteriomente,establecemos como datos faltantes el periodo que deseamos predecir en este caso 2019 .
db <- db %>%
mutate(n=ifelse(annus == 2019 & prov=="LIMA" ,NA,n))
db.n <- db.lima %>%
mutate(n=ifelse(annus == 2019,NA,n)) Posteriormente, de nueva cuenta empleamos la funcion inla() con todos los parametros complementarios que detallamos. Sin embargo, un punto importante a considerar es que el proceso que modelamos es un proceso de conteo por lo tanto el supuesto respecto de la distribucion que emplearemos sera el de una binomial negativa family ="nbinomial" . Cabe mencionar que podemos ver la lista completa de distribuciones posibles para la variable de respuesta empleando :inla.list.models("likelihood")
# Asumiendo una Distribucion apriori para los años :"ar1"
lima.m.m1 <- inla(n ~ 1 + temperature + f(annus,model = "ar1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
verbose = F,
family = "nbinomial",
data = db.n
)
# Asumiendo una Distribucion apriori para las semanas :"rw1"
lima.m.m2 <- inla(n ~ 1 + temperature+ f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n
)
#
lima.m.m3 <- inla(n ~ 1 + temperature + f(annus,model = "ar1") + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n
)Una vez ajustados los modelos podemos recolectar informacion relevante de su ajuste. En este sentido el objeto generado por la función inla() es una lista que contiene todos los resultados los cuales se corresponden con las opciones(parametros) utilizadas durante la ejecucion de inla(). En este caso recurrimos a objetoInla$summary.fitted.values para recolectar informacion de los valores ajustados.
fit.m.m1 <- lima.m.m1$summary.fitted.values$mean
fit.m.m2 <- lima.m.m2$summary.fitted.values$mean
fit.m.m3 <- lima.m.m3$summary.fitted.values$mean
n.m <- db.lima$nAsimismo procemos a asignar toda la informacion de los valores ajustados de los modelos a un objeto lista. Posteriormente empleamos melt para crear una estructura de base de datos del tipo larga adecuada para graficar .Asi como recolectamos informacion de los intervalos de credibilidad de los valores ajustados empleando : lima.m.m1$summary.fitted.values$0.025quant y lima.m.m1$summary.fitted.values$0.975quant,
datos <- list("AR(1)"=fit.m.m1,"RW(1)"=fit.m.m2,"AR(1)RW(1)"=fit.m.m3) %>%
as.data.frame() %>%
melt( variable.name = "modelo",
value.name = "fit") %>%
group_by(modelo) %>%
mutate(actual = n.m,
date = db.n$date) %>%
ungroup() %>%
mutate(min_inter = c(lima.m.m1$summary.fitted.values$`0.025quant`,
lima.m.m2$summary.fitted.values$`0.025quant`,
lima.m.m3$summary.fitted.values$`0.025quant`),
max_inter = c(lima.m.m1$summary.fitted.values$`0.975quant`,
lima.m.m2$summary.fitted.values$`0.975quant`,
lima.m.m3$summary.fitted.values$`0.975quant`)
)Posteriormente, graficamos los datos asi como el ajuste de los modelos propuestos con su respectivo intervarlo de credibilidad empleando la funcion ggplot() del paquete ggplot2
datos %>%
ggplot() +
geom_line(aes(x=date,y=actual),lwd=1.5) +
geom_line(aes(x=date,y=fit),lwd=1.5,color="#011f4b",linetype = "dashed") +
geom_ribbon(aes(x=date,y=fit,ymin=min_inter, ymax=max_inter),
alpha=0.2, #transparency
fill="#011f4b",
color="#011f4b",lwd=1.5) +
facet_wrap(vars(modelo)) +
theme_bw()+
xlab(" ") +
ylab("Numero de casos")+
geom_vline(xintercept = as.Date("2019-01-07"), linetype="dotted",
color = "black", size=3)+
theme(plot.title = element_text(hjust = 0.5,size = "20pts"),
strip.background =element_rect(fill="#011f4b"),
strip.text = element_text(color = "white",face = "bold",size = "35")) 
2.2.2 Performance
Para el calculo de metricas en la búsqueda de evaluar los modelo propuestos emplearemos el paquete yardstick . En particular empleamos la funcion metric_set para seleccionar la metricas que deseamos calcular (mae,rmse,msd). Agrupamos los datos por cada modelo para realizar el calculo a ese nivel, filtramos por el año de analisis de la prediccion 2019 y finalmente empleamos la funcion perform.metrics, para lo cual detallamos los parametros truth y estimate que recogen los datos verdaderos y ajustados, respectivamente.
#Metricas a emplear
perform.metrics <- metric_set(mae, rmse,msd)
#Calculo de las metricas en el año de interes de la prediccion
tbl.yrd.m <- datos %>%
filter(date>=as.Date("2019-01-07")) %>%
group_by(modelo) %>%
perform.metrics(truth = actual,
estimate = fit)Finalmente, el objeto resultado de dichas métricas tiene estructura del tipo larga por lo que para efectos de reportar estos resultados en una tabla. En primer lugar la estructura data.frame es convertida a ancha empleando la funcion pivot_wider la cual recolecta como nuevos nombres de columnas(names_from) los niveles de la variable .metric y estas a su vez se rellenan con los valores(values_from) de .estimate. Finalmente empleamos del paquete gt la funcion gt() para reportar una tabla personalizada .
#tabla de resultados
tbl.yrd.m %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de Precisión fuera de la muestra"),
subtitle ="out-of-sample accuracy metrics") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de Precisión fuera de la muestra | |||
|---|---|---|---|
| out-of-sample accuracy metrics | |||
| modelo | mae | rmse | msd |
| AR.1. | 89.32137 | 105.10674 | 88.29872 |
| RW.1. | 74.96485 | 91.60271 | 74.23950 |
| AR.1.RW.1. | 89.23152 | 101.74050 | 89.23152 |
De las estructuras temporales el mejor modelo en terminos de prediccion fuera de la muestra (forecasting out-of-sample) fue el modelo que asume un paseo aleatorio de orden 1 para dicha estructura `rw1`` .
2.3 Prediccion: Enfoque Espacio-Temporal
El enfoque espacio-temporal, implica complementar el enfoque anterior complejizando el modelo tomando en cuenta la dimesion espacial. Para ello requerimos incluir efectos aleatorios espaciales para controlar por tal dimension de analisis. Cabe mencionar que los modelos espaciales areales empleando INLA, requieren necesariamente la creacion de un grafo(graph) o una estructura de datos espaciales que capture la disposicion de las areas como una matriz de pesos espaciales.
\[\begin{align*} y_{t} = S_{t}+T_{t}+\nu_{t}+u_{t}+e_{t} \end{align*}\]
donde:
\(\nu_{t}\) : efectos estructurados
\(u_{t}\) : efectos no estructurados
Los 2 nuevos componentes \(\nu_{t},u_{t}\) pretenden tomar en consideracion en la prediccion la correlacion espacial presente en los datos . Al respecto :
Los efectos no estructurados corresponden a efectos aleatorios que pretenden controlar caracteristicas no observadas de cada area bajo estudio. Para modelarlas podemos emplear
f(area, model = "iid")Los efectos estructurados son efectos aleatorios que toman en cuenta explicitamente la estructura espacial
- Pueden ser modelados en INLA de diversas maneras:
- Empleando efectos espaciales besag :
f(area, model = "besag")
- Empleando efectos espaciales propios de besag :
f(area, model = "besagproper")
- Empleando efectos espaciales besag :
- Pueden ser modelados en INLA de diversas maneras:
2.3.1 Preprocesamiento
Dado lo expuesto anteriormente, recolectamos la informacion geografica de Peru expresada en un shapefile, a nivel provincial . Para lo cual ademas de llamar al shapefile, al agrupar los datos a nivel provincial y emplear la funcion summarize buscando condensar los poligonos del shapefile al nivel deseado
data("Peru")
peru.prov<-Peru %>%
group_by(reg,prov) %>%
summarise()Posteriomente, creamos una estructura de vecindad a partir del shapefile , y a su vez de dicha estructura de vecindad creamos la matriz de pesos espaciales
# Creando la estructura de vecinos mas cercanos
peru.adj <- poly2nb(peru.prov)
# Pesos espaciales
W.peru <- nb2mat(peru.adj, style = "W") Finalmente, realizamos la limpieza de provincias inexistentes; asi como la creacion de un id por cada provincia, para asi poder identificar cada poligono. Esto ultimo representa un requisito extra para ajustar modelos espaciales en INLA
db.m.sp<- db %>%
group_by(prov) %>%
filter(!prov %in% c("EXTRANJERO","ARICA")) %>%
mutate(id.sp=cur_group_id())2.3.2 Ajuste y prediccion
Ajustamos el modelo, siguiendo los pasos antes indicados
peru.m.m5 <- inla(n ~ 1 + f(annus,model = "linear") +
f(week,model="rw1") +
f(id.sp, model = "besag",
graph=W.peru)+
temperature
#f(inla.group(temperature),model="rw1")
#temperature
,
data = db.m.sp,
control.compute = list(dic = TRUE,
waic = TRUE,
cpo = TRUE),
control.predictor = list(link = 1)
)De nueva cuenta recolectamos la informacion del objeto resultante de emplear inla()
db.m.sp.fit <- db.m.sp %>%
ungroup() %>%
mutate(fit =peru.m.m5$summary.fitted.values$mean,
min_inter=peru.m.m5$summary.fitted.values$`0.025quant`,
max_inter=peru.m.m5$summary.fitted.values$`0.975quant`)
db.lima.m.fit <- db.m.sp.fit %>%
filter(prov=="LIMA")
fit.m.m1 <- lima.m.m1$summary.fitted.values$mean
fit.m.m2 <- lima.m.m2$summary.fitted.values$mean
fit.m.m3 <- lima.m.m3$summary.fitted.values$mean
fit.m.m4 <- db.lima.m.fit$fit
min.m.m4 <- db.lima.m.fit$min_inter
max.m.m4 <- db.lima.m.fit$max_inter
n.m <- db.lima$n
datos.m <- list("AR(1)"=fit.m.m1,"RW(1)"=fit.m.m2,
"AR(1)RW(1)"=fit.m.m3,"AR(1)RW(1) with Spatial effects"=fit.m.m4) %>%
as.data.frame() %>%
melt( variable.name = "modelo",
value.name = "fit") %>%
group_by(modelo) %>%
mutate(actual = n.m,
date = db.n$date) %>%
ungroup() %>%
mutate(min_inter = c(lima.m.m1$summary.fitted.values$`0.025quant`,
lima.m.m2$summary.fitted.values$`0.025quant`,
lima.m.m3$summary.fitted.values$`0.025quant`,
min.m.m4),
max_inter = c(lima.m.m1$summary.fitted.values$`0.975quant`,
lima.m.m2$summary.fitted.values$`0.975quant`,
lima.m.m3$summary.fitted.values$`0.975quant`,
max.m.m4)
) y de nueva cuenta graficamos los resultados del modelo para comparar su ajuste respecto de los datos verdaderos .
2.3.3 Performance
Finalmente, seguimos los mismos pasos que en la sección del enfoque temporal.
#Metricas a emplear
perform.metrics <- metric_set(mae, rmse,msd)
#Calculo de las metricas en el año de interes de la prediccion
tbl.yrd.m <- datos.m %>%
filter(date>=as.Date("2019-01-07")) %>%
group_by(modelo) %>%
perform.metrics(truth = actual,
estimate = fit)
#tabla de resultados
tbl.yrd.m %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de Precisión fuera de la muestra"),
subtitle ="out-of-sample accuracy metrics"
) %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de Precisión fuera de la muestra | |||
|---|---|---|---|
| out-of-sample accuracy metrics | |||
| modelo | mae | rmse | msd |
| AR.1. | 89.32137 | 105.10674 | 88.29872 |
| RW.1. | 74.96485 | 91.60271 | 74.23950 |
| AR.1.RW.1. | 89.23152 | 101.74050 | 89.23152 |
| AR.1.RW.1..with.Spatial.effects | 76.11171 | 87.56624 | 72.87075 |
El modelo espacial como se puede ver si bien no recoge los patrones temporales del todo bien, en promedio en terminos de precision es el mejor modelo resultante pues no solo obtiene predicciones con intervalos de credibilidad más acotados sino tambien supera a los otros modelos en terminos de la prediccion fuera de la muestra(out-of-sample). Finalmente cabe mencionar que estos resultados pueden ser ampliamente mejorados con la especificacion adecuada del predictor linear usando estructuras apriori más adecuadas a los datos o quiza otras más complejas;asi como considerando otras variables de control.